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ABSTRACT 

Context. Radiative shocks play a dominant role in star formation. The accretion shocks on the first and second Larson's cores involve 
radiative processes and are thus characteristic of radiative shocks. 

Aims. In this study, we explore the formation of the first Larson's core and characterize the radiative and dynamical properties of the 
accretion shock, using both analytical and numerical approaches. 

Methods. We develop both numerical radiation-hydrodynamics calculations and a semi-analytical model that characterize radiative 
shocks in various physical conditions, for radiating or barotropic fluids. Then, we perform ID spherical collapse calculations of the 
first Larson's core, using a grey approximation for the opacity of the material. We consider three different models for radiative transfer, 
namely: the barotropic approximation, the flux limited diffusion approximation and the more complete Ml model. We investigate the 
characteristic properties of the collapse and of the first core formation. Comparison between the numerical results and our semi- 
analytical model for radiative shocks shows that this latter reproduces quite well the core properties obtained with the numerical 
calculations. 

Results. The accretion shock on the first Larson core is found to be supercritical, i.e. the post and pre-shock temperatures are equal, 
implying that all the accretion shock energy on the core is radiated away. The shock properties are well described by the semi- 
analytical model. The flux limited diffusion approximation is found to agree quite well with the results based on the Ml model of 
radiative transfer, and is thus appropriate to study the star formation process. In contrast, the barotropic approximation does not 
correctly describe the thermal properties of the gas during the collapse. 

Conclusions. We have explored and characterized the properties of radiative shocks typical of the formation of the first Larson's 
core during protostellar collapse, using both radiation-hydrodynamics numerical simulations and a semi-analytical model, except, 
for this latter, in the case of subcritical shocks in an optically thin medium. We show that a consistent treatment of radiation and 
hydrodynamics is mandatoiy to correctly handle the cooling of the gas during the core formation and thus to obtain the correct 
mechanical and thermal properties for this latter We show that the flux limited diffusion approximation is appropriate to perform star 
formation calculations and thus allows a tractable and relatively correct treatment of radiative transfer in multidimensional radiation- 
hydrodynamics calculations. 

Key words. Stars: formation - Methods : analytic, numerical - Hydrodynamics - Radiative transfer 

1. Introduction acteristic of the first core initial energy content. This entropy 

content is kept roughly constant during the nearly adiabatic sub- 

Star formation involves a large variety of complex physical ^^^^^^^ ^^^^^^ ^^jj^p^^ formation of the second 

processes. Among them, radiative transfer appears to be one of ^^^.^ accretion shock on the first core is thus of prime im- 

the most important ones, since it governs the behaviour of the portance, since it determines the entropy level of the following 

H eariiest phases of the collapse and the formation of the so-called ^^^^^^ formation, up to the formation of the protostar it- 

. . first core (ILarsonlll969 l). It is well established that the molec- ^^^^ 
ular gas experiences different thermal regimes, from isother- 
mal to adiabatic, during the earliest phase of the collapse, via 

the coupling between gas, dust and radiation (e.g. Larson l l 19691: A radiative shock is a shock that is so strong that it emits 

iTscharnuter & Winkledll979l iMasunaea et al.lll998l) . radiation which in turn affects the hydrodynamic behaviour of 

Since the pioneering works of Winkler & Newman and the flow. The equations of radiation-hydrodynamics (RHD) must 

Tscharnuter (Winkler & Newman 1980; Tscharnuter 1987), rel- thus be solved consistendy to properly explore this kind of a 

atively few studies have been devoted to the accretion shock on process. Radiative shocks are common in astrophysics, e.g. in 

the first Larson's core. This shock is a radiative shock, i.e. radi- star formation regions or around supernovae remnants. They 

ation can escape from the infalling material. Once this infalling are well studied in the literature, from th e theoretical, experi- 

gas has been shocked and has reached an optical depth of 1, it mental and numerical points of view (e.g. iBouquet et al.ll2000l : 

becomes more or less adiabatic, at a given entropy level char- JDrake 2007; Gonzalez et al. 2009). The large vari ety of radia- 
tive shocks makes it difficult to try to classify them (lDrakel2005l; 
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In this work, we study radiative shocks from both numeri- 
cally and analytically. We focus on the properties of the accre- 
tion shock during the formation of the first prestellar core. In 
the first part of the paper, we recall the main properties of ra- 
diative shocks. Then, considering the jump relations for a radiat- 
ing fluid, we focus on some particular kinds of radiative shocks, 
with upstream and downstream material characterized by dif- 
ferent optical depth properties. We also consider the case of a 
shock taking place in a barotropic material, since the barotropic 
approximation is the simplest way to charac terize the thermal 
evolution of the gas during the collapse (e.g. ICommer^on et al.l 
[2008 ). 

In the second part of this work, we study the impact of var- 
ious models of radiative transfer on the protostellar collapse. 
Radiation hydrodynamics plays here a crucial role, for instance 
by evacuating the compressional energy, leading to a nearly 
isothermal free-fall collapse phase. Radiation transfer also has 
a dramatic impact on the accretion shock, with the infalling gas 
kinetic energy being either converted into internal energy in the 
static adiabatic core or radiated away. This ratio between the en- 
ergy accreted onto the star and the one radiated away is a key 
quantity, as it ultimately determines the entropy content of the 
forming star. Using a ID spherical code, we compare results of 
calculations done with a barotropic EOS, a flux limited diff'usion 
approximation, and a Ml model for radiative transfer, assum- 
ing grey opacity for the infalling material. ID calculations retain 
the important virtue of allowing detailed and complex physical 
processes to be considered in dense core collapse calculations, 
which is not the case in a multidimensional approach. In conse- 
quence, a ID approach enables us to characterize the impact of 
these various physical processes on the results, allowing a better 
characterizat ion of the validity of s i mplifi ed multidimensional 
studies (e.g. Commer9on et al .1120101 1201 ll) . 

This paper is organized as follows. In the first section, we 
recall the basic properties of radiative shocks and develop semi- 
analytical models that can be applied in some particular cases 
for a radiating fluid and a barotropic fluid. In Sect. 3, we detail 
the numerical method and the physics inputs used in our calcu- 
lations. In Sect, m we derive the first Larson's core properties 
from the numerical calculations, using various approximations 
for the radia tive transfer, and compa re the results with the ones 
obtained by iMasunaga et al.1 (Il998h . We also present an origi- 
nal semi-analytical model that allows a simple interpretation of 
the numerical results. Eventually, Sect. 5 summarizes our main 
results. 



2. Radiative shock - A semi-analytic model 

2.1. A qualitative picture of radiative sliocks 

Depending on the shock's strength, radiative shocks be- 
long to two groups: the subcritical shocks and the supercritical 
shocks. As the strength of a shock increases, the postshock tem- 
perature T2 rises, producing a radiative flux of order crT^ that in- 
creases very rapidly. This flux penetrates the upstream material 
and preheats this latter to a temperature T immediately ahead 
of the shock front (radiative precursor) that is proportional to the 
incident flux. T increases rapidly with the shock strength and 
eventually becomes equal to 7^2. A shock with T < T2 is called 
a subcritical shock. Because the material entering the shock is 
preheated, the postshock temperature T+ exceeds its asymptotic 
equilibrium value T2, and decays downstream as the material 
cools by emitting photons that propagate across the shock (see 
Fig. 1 for various shock characteristics). 



For stringer shocks, the preheating becomes so important 
that the preshock temperature equals the postshock equilibrium 
temperature T2. The shock velocity at which the postshock and 
preshock temperatures are equal, defines the critical shock. For 
higher shock velocities, T cannot exceed T2, the excess energy 
forces the radiative precursor further into the upstream region 
with a temperature close to T2. Pre- and post-shock temperatures 
are equal, the supercritical shock is thus isothermal. Radiation 
and matter are still out of equilibrium in some part of the precur- 
sor, but come to equilibrium when temperature approaches T2. 
The upstream kinetic energy is radiated at the shock. 

Most of the early work on radiative shocks has be en de- 
scribe d in lZel'Dovich & Raizerl (il967) and Mihalas & MihalasI 
where readers can find the basic equations of the front 
structure, radiative precursor extension, etc... 

2.2. Jump relations for a radiating material 

Consider th e jump relations (Rankine Hugoniot) for a radi- 
ating flow (see Mihalas & Mihalasl ll984t IZel Dovich & Raizerl 
1967; Drake 2007.) . Conservation of mass, momentum and en- 
ergy yield 

piM] — P2U2 = m, (1) 
piM^ + Pi+ = p2ul +P2+ P,2, (2) 

in + piufj +Frl +Ui (Eyi + f ri) = 

m {h2 + P2ul) + F,2 + U2 (E,2 + Prl) , (3) 

where subscripts "1" and "2" denote respectively the upstream 
and downstream states, all radiation quantities are estimated in 
the comoving frame, h corresponds to the gas specific enthalpy 
(h - e + Pip, with e - pe - P I iy - \)) and m is the mass flux 
through the shock. In comparison with the hydrodynamical case, 
the pressure is now the total pressure, i.e. the gas plus radiation 
contributions, P+Pi, and the specific enthalpy is the total specific 
enthalpy, e + (P + P^ + Ei)/p. The radiation energy Ei and the 
radiation pressure P^ are important only at high temperatures or 
low densities, whereas the radiation flux, F^, plays a fundamental 
role in all radiative shocks. 

Contrarily to the hydrodynamical case, the system of equa- 
tions ([T]i, (O and (O can not be solved explicitly. We need to 
make some assumptions on both the upstream and downstream 
materials. In what follows, we distinguish two cases: i) the up- 
stream and downstream materials are opaque, ii) the upstream 
material is optically thin and the downstream material remains 
opaque. 

2.2.1 . Radiative sliock in an optically thick medium 

This is the most studied case (e.g. lMihalas&Mihalaslll984l: 
lDrakell2007h . It occurs, for instance, in shocks in stellar interi- 
ors, where matter is both dense and hot. At a sufficiently large 
distance from the front, matter and radiation are in equilibrium 
and both the downstream and upstream materials are opaque 
[Fi-i = Fi2 = 0). Any radiation crossing the front from the 
hot downstream material into the cool upstream material is re- 
absorbed in the radiative precursor, into which it propagates by 
diffusion. Outside this diffusion layer, the flux vanishes and we 
have 

til + piufj+ui (£1-1 + P,i) = m {h2 + P2ul)+U2 (E,2 + P12) -(4) 

Since matter and radiation are in equilibrium and the mate- 
rial is opaque, we have E^ = 3Pi = aj^T'^. Defining the compres- 
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Fig. 1. Temperature and d ensity profiles in a subcritical shock (left) and a supercritical shock (right). Adapted from 
IZeF Dovich & Raized (Il967h . 



sion ratio r = pi/pi - mi/m2, we can rewrite equations ^ and 
(O in the non-dimensional form 



and 



2 1'" ^i^(n-i) + Q'i 



2-1 



y-l\r / \ 



(5) 



(6) 



where n = Pa/^"!, = (l/3)flRr^/Pi and Mi is the hydro- 
dynamic Mach number, Aii = Ui/cs\ = Mi/(yPi/pi)'^2 xhe 
coupled equations (|5]l and (|6]l are solved numerically to get the 
variations of H and r as function of Mi and . 

The compression ratio rises from r = 4 for a pure hydrody- 
namic shock with y = 5/3tor = 7, which corresponds to the 
limiting compression ratio for a gas with y = 4/3, i.e pure radia- 
tion, a gas made of photons. The stronger the shock, the greater 
the radiative effects, even for a very small initial ratio of radiative 
pressure to gas pressure. When the upstream radiative pressure 
increases, the shock becomes immediately radiative since radia- 
tive pressure increases as T'^, whereas the gas pressure increases 
only li nearly with temperature. As shown in Mihalas & Mihalas 
(|1984|) . for a strong radiating shock, since the compression ratio 
is fixed, the temperature ratio increases only as M\^^, whereas it and 
rises as M^ in a non-radiating shock. Note that if Pi/P < 10"^, 
the non-radiating fluid approximation becomes valid for an up- T2 
stream Mach number All < 10. ~ 



2.2.2. Radiative sliock witli an optically thin upstream 
material 

Suppose now that the shock is propagating in an optically 
thin material with an opaque downstream region, as in the case 
of star formation (e.g., Calvet & Gullbring 1998) . In the low 
mass star formation context, the radiative energy and pressure 
can be neglected compared to the gas internal energy for the first 
core accretion shock. 

Let us consider the material and radiative quantities at the 
discontinuity, outside the spike in gas temperature (i.e., between 
regions with subscripts "2" and "-"). We have 



P2U2 = m. 



p-U_ + — P2U2 + P2, 
til (h- + p-u^^ = m {h2 + P2M2) + AFr, 



(7) 
(8) 

(9) 



where AFi = F^2 - Fr-- We consider the case of a non zero net 
flux \Fi2 - Fr-\ across the shock. The jump relations, derived 
from the conservation equations ([T), (|2]l and (|3]l, read 



P2 

P- 



[(y + i)f 2 + (y - i)P-] U2 + 2(y - \)^F, 

[(y + 1)P_ + (y - Y)P2\ U2 



[(y + 1)P_ + (y - Y)P2\ M- - 2(y - 1)AF, 
[(y ^ \)P2 + (y - 1)P-] M- 



(10) 



(11) 
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These two relations show that radiative energy transport (i.e. the 
radiative flux) across a shock can significantly alter the density, 
temperature and velocity profiles of the flow. Both upstream and 
downstream materials are affected over distances which depend 
on the material opacity. As mentioned before, the stracture of a 
radiative shock is as follows: the upstream material is preheated 
by a radiation precursor while the downstream material is cooled 
by radiative losses. 

As for the previous opaque case for which we use an iterative 
solver, the downstream quantities can not be derived analytically 
from the conservation relations for given upstream conditions. 
The radiative flux has to be known and the result depends on the 
upstream flow. 

The simplest case to study is a supercritical shock, where 
r_ = T2 and then /i_ = /z2. Using the same adimensional param- 
eters as for the opaque case, we have 



and 



1 



= n- 1, 



y-\\r I P^u^ 



Note that since the shock is isothermal, r -W and 



(12) 



(13) 



(14) 



Eventually, the radiative flux discontinuity, normalized to the up- 
stream kinetic energy, is given by: 



AF, _ y^At - 1 
Q.Spiul ~ y^Mt ' 



(15) 



where X thus represents the amount of incident kinetic energy 
radiated away at the shock front. 

Figure |2] shows the evolution of the adimensional compres- 
sion ratio, r, and flux discontinuity, X, as a function of the up- 
stream Mach number for a supercritical shock. As seen in the 
figure, the compression ratio does not saturate, in contrast to the 
case of an opaque material. As shown by equations ( fTOl l and ( fTTT i. 
a non-zero flux across the front, with F,2 > F,-, increases the 
density jump whereas it decreases the temperature jump. In or- 
der to have an isothermal shock at low Mach number, Al < 2, 
the downstream velocity is determinant since in that case the 
upstream kinetic energy is not entirely radiated away. In proto- 
stellar collapse calculations, the characteristic Mach number at 
the first core accretion shock is ~ 3, implying that almost all 
the kinetic energy (> 99%) is radiated away at this stage. 

Including radiative energy and radiative pressure in equa- 
tions ^ and (|9]l does not affect the dependence of the compres- 
sion ratio upon the upstream Mach number. On the other hand, 
the dependence of the amount of kinetic energy radiated away 
becomes stronger: 



X 



AF, y^Af - 1 



0.5p-ul y^At 



(1 -H4a'i). 



(16) 



The same analysis cannot be carried out for a subcritical 
shock with an optically thin upstream medium, since, in that 
case, no constraint allows us to close the system of equations 
(O and (|9]l. One needs a prescription on the radiation flux or 
luminosity in the upstream material to close the system. 
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Fig. 2. Compression ratio r and fraction of kinetic energy radi- 
ated away at the shock, X, as a function of the upstream Mach 
number for a supercritical shock, in the case of an opticaUy thin 
upstream material. 

2.3. Super- or sub-critical shock ? 

In this section, we characterize in which regime the radiative 
shock takes place, depending on the upstream and downstream 
material properties. 

2.3.1. Opaque material 

In an opaque material, the simplest case is the one of a super- 
critical shock, in which the preshock gas ahead of the disconti- 
nuity is heated up to the postshock temperature (see §2. 1 .2). The 
radiative energy absorbed in the upstream region is used only to 
raise the gas temperature (Zel'Dovich & Raizer 1967, p. 536). 
Assuming, for sake of simplicity, that the gas is neither com- 
pressed nor slowed down, that the radiative pressure is negligible 
compared to the the gas pressure (valid for early low mass star 
formation stages), and that the upstream internal energy is negli- 
gible far from the shock (region with subscript "1"), the equation 
of energy reads 



F, = u-p-e{TS), 



(17) 



with e(T) = e - If we evaluate the flux iust at the 

discontinuity, we find the maximum preheating temperature T_ 
from 

\F,\ ~ crr^ = M_p_e(r_). (18) 

We can now determine the critical temperature at which 
T equals T2 

o-Ti = M_p_e(r„), (19) 
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so that 



are simply 



1/3 



(20) 



(y- l)fimHO-) 
which defines the supercritical shock condition. 

2.3.2. Optically thin upstream material 

In the case of an optically thin upstream material, the afored- 
erived criterion for a supercritical shock is simply derived by as- 
suming that all the upstream kinetic energy is radiated away at 
the shock, which yields 

7^0.= • (21) 



2.4. Estimate of the preshock temperature 

In this section, we calculate the preshock temperature as a 
function of the upstream quantities, for various natures of the 
shock. 



2.4.1 . Supercritical shock with an optically thin or thick 
upstream material 

This is the simplest case since the results are independent of 
the optical depth of the upstream material. Once the upstream 
density is known, the velocity and the temperature are set. When 
the upstream gas is not compressed (ui = m_) and for a strong 
shock (Ai > 2 so that X ~ 1), it is easy to get the shock tem- 
perature Ts - - T2 since all the upstream kinetic energy is 
radiated away: 



(O.Spiu 



3\ 



1/4 



(22) 



2.4.2. Subcritical shock with an optically thick upstream 
material 

In that case, the shock temperature is fixed by the up- 
stream velocity. The equ ilibrium postshock temperature is ( 
iMihalas & MihalasI (11984 -): 



T2 = 



2(7- l)ut 



(23) 



Kiy + 1)2 ' 

with = I until. The preshock temperature T is estimated as 



r-1 



(24) 



which indicates that at any location ahead of the shock, all the 
radiative energy going across this location is absorbed and heats 
up the gas. 

2.5. Jump relations for a barotropic gas 

A widely used approximation to handle radiative transfer 
during the first stages of star formation is the barotropic approxi- 
mation (e.g. Commercon et al. 2008). For a barotropic material, 
the temperature and pressure depend only on the density. The to- 
tal energy is thus not conserved. In this case, the jump relations 



piui - P2U2 = m, 
p\u\ -hPi = P2u\ +P2, 



(25) 
(26) 



with the barotropic equation of state (EOS) as a closure relation 



P cc p 



1 + 



Pad 



y-1 



(27) 



where is the critical density at which the gas becomes adia- 
batic (see Sect. I3.4.T] ). Using the same adimensional parameters 
as previously, we get 



21-1. = n-i. 



Using the barotropic EOS yields 



n 



Pi P2 1 + 



,(y-i) 



1 + /r-i)^(r-i) 



Pi 1 + r 



,(r-i) 
I 



1 + r 



.(r-i) 



(28) 



(29) 



where ri = pi/pad and r2 = pi/Pad- Eventually, we get r by 
solving 



1 +r' 



,(y-i) 



,r-l 



= (r- l) + r]'\r^ - 1) 



(30) 



One can also derive an equivalent luminosity, assuming a 
shock in a perfect gas with the same jump properties. From equa- 
tion (|9]l, we have 



X = 



AFr 



1 



O.Spiu] 7-1 



(31) 



where r and 11 are given by (l30l l. 

Figure [3] illustrates the evolution of the adimensional quan- 
tities r and X as a function of the upstream Mach number for a 
barotropic material, for density ratio ri = 0.1, 1, and 10, which 
represent transition values between the optically thin (isother- 
mal) and thick regimes (adiabatic). For comparison, we also plot 
again the evolution of r and X for a supercritical shock with an 
optically thin upstream material. The compression ratio is not 
bound by a limiting value as the upstream Mach number in- 
creases, as for the case of a supercritical shock with an optically 
thin upstream material, but in contrast to the case of an opaque 
material. The compression ratio decreases as ri increases (note 
that for ri = 10, the material should be totally optically thick in 
the context of star formation). The shock is never found to be 
supercritical (X ~ 1), even though almost all the incident kinetic 
energy is assumed to be radiated away at high upstream Mach 
number. This clearly shows that a barotropic EOS approxima- 
tion cannot treat properly radiative shocks. Compared to the case 
of a supercritical shock with an optically thin upstream which is 
relevant for the accretion shock on the first Larson core in the 
transition regions between the isothermal and adiabatic regime, 
the compression ratio and the amount of energy radiated away 
are underestimated. 

The specific entropy jump at the shock is expressed as 



As = S2 - = Cy 



(32) 



where s is the specific entropy and Cy = ks/(jimii(y - 1)). 
The entropy jump for a barotropic fluid is thus proportional 
to n/r''. Figure |4] shows the evolution of log(n/r'') for a 
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Fig. 3. Evolution of the compression ratio, r, and of the amount 
of kinetic energy radiated away, X, as a function of the upstream 
Mach number for a barotropic material, with density ratio ri - 
0.1 (black), 1 (red), and 10 (blue). The dotted lines represents the 
evolution for a supercritical shock with a optically thin upstream 
material (c.f. Fig.|2]l. 

barotropic fluid (with r\ - 0.1, 1, and 10) and for a shock with 
an optically thin upstream material. The entropy jump at the 
shock is overestimated with a barotropic law. This would lead 
to an incorrect initial entropy level and profile for new born 
protostars. The entropy is indeed a key quantity for pre-main 
sequence evolutionary models as it entirely determines the 
mass-radius relationship of an adiabatic object, like the Larson 
cores. 

In this section, we have shown that the downstream quan- 
tities of radiative shocks strongly depend on the nature of the 
shock and on the model used to describe radiation transport. In 
the following sections, we study in detail the case of a particular 
radiative shock: the accretion shock on the first Larson core, in 
the context of star formation. 

3. 1D spherical numerical calculations - Method. 

3.1. introduction and previous work 

In this subsection, we introduce our numerical method and 
basics concepts, for a good understanding of the following part 
of this work. 

The two companion papers Ma sunaga et al.l (11998 '). and 
iMasunaga & I nutsuka (2000) present a very extensive study of 
ID protostellar collapse. The first paper focusses on the for- 
mation and the properties of the first core, using a frequency 




S -1.0 - ■ 

M 
O 

-1.5- - 

-2.0 [ 

1 10 
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Fig. 4. Evolution of the entropy jump, Ai oc II/ r^, as a function 
of the the upstream Mach number for a barotropic material, with 
density ratio r\ - 0.1 (black), 1 (red), and 10 (blue). The dotted 
line represents the evolution for a supercritical shock with an 
optically thin upstream material. 

dependent RHD model, while the second paper addresses the 
formation of the second core These authors use a ID spherical 
code in the Lagrangean frame. Moment equations of radiation 
in the comoving frame are solved following IStone et al.l (Il992l) . 
iMasunaga et alJ (Il998h use a variable Eddington tensor factor 
(VETF) that retains a frequency dependence, whereas grey opac- 
ities are used to calculate the coupling between matter and radi- 
ation and the work of t h e radi ative pressure. 

In IMasunaga et alj ( Il998h . the cores initially have uniform 
density and temperature, and a radius adjusted so that the 
cloud should be initially slightly more massive than the Jeans 
mass. Initial masses range from 0.1 M© to 3 M©. They find 
that whatever the initial cloud mass, the first core radius Rf^ is 
almost constant, with Rf^ ~ 5 AU, and the first core mass is 
Mfc ~ 0.05 Mq. In their study, Rf^ is defined as the point where 
the gas pressure is balanced by the ram pressure of the infalling 
envelope. 

To compare our numerical results with analytical theories, it 
is useful to define some measurable quantities. A useful one is 
the mass accretion rate M. For a ID spherical model, the mass 
accretion rate is simply M = An r^pu. In the theory, the accretion 
rate is generally defined as (e.g. lStahler et al.l[l980h 

M = (33) 

where ^ is a dimensionless coefficient and Cso the isothermal 
sound speed. ^ is estimated assuming that a non-magnetic, non- 
rotating cloud, whose initial state corresponds approximately to 
a balance between thermal support and self-gravity, has com- 
parable free-fall time, % ~ R^^^ /(GMc)^^^, and sound-crossing 
time, fsc ~ Rr/C s (). The mass accretion rate is thus given by 
M ~ MJtff. ' Shul (11977b obtained ^ = 0.975 for the expansion- 
wave solution , whereas the dynami c Larson-Penston solution 
yields ^ ~ 50 (lLarsonlll969HPenston|[l 969). 

Another quantity directly comparable to the theory is the ac- 
cretion luminosity, usually expressed as 

z.r = 5^. (34) 

where Mfc corresponds to the mass of the accreting core. This 
accretion luminosity thus corresponds to the case where all the 
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infalling kinetic energy is radiated away, i.e. to a supercritical 
radiative shock. The accretion luminosity can be directly mea- 
sured in numerical calculations, from the mass accretion rate, 
and compared with the intrinsic effective luminosity emerging 
from the first core. 



to an adiabatic, monoatomic gas with adiabatic exponent y - 
5/3. Molecular hydrogen behaves like a monoatomic gas until 
the temperature reaches several hundred Kelvin since the rota- 
tional degrees of freedom are not excited at lower tem peratures 
(IWhitworth & Clarke|[l997t iMasunaga & Inutsukal200Q[) . 



3.2. Numerical method 



3.4.2. Moment Models 



We use a ID full Lagrangean version of the co de devel- 
oped by Chieze and Audit (see ' Audit et al] 1200 2*). that in- 
tegrates the equations of grey radiation hydrodynamics under 
three different assumptions (in order of increasing complexity 
order): a barotropic EOS, the flux limited diff'usion approxima- 
tion (FLD) and the Ml model. The RHD equations are integrated 
in their non-conservative forms using finite volumes and an ar- 
tificial viscosity scheme in tensorial form. The RHD equations 
are integrated with an implicit scheme in time, using a standard 
Raphson-Newton iterative method. 

3.3. General assumptions for the radiation field 

The first approximation in our calculations is to consider 
a grey radiation field, i.e. only one group of photons, inte- 
grated over all frequencies. We also neglect scattering and as- 
sume Local Thermodynamical Equilibrium (LTE) everywhere. 
The second assumption concerns the coupling between radia- 
tion and hydrodynamics, for which we write the RHD equations 
in the comoving frame. 

3.4. Models for the radiative transfer 

According to p revious studie s using accurate model for the 
radiation field (e.g. lMasunaga et al. 1998), it is well established 
that the molecular gas follows two thermal regimes during the 
first collapse. At low density, the gas is able to radiate freely and 
to couple with the dust. This is the isothermal regime, where the 
Jeans mass decreases with density. Then, when the gas becomes 
denser, typically p > ~ 1 x lO"'-' g cm"-', the radiation is 
trapped within the gas that begins to heat up adiabatically. The 
Jeans mass then increases with density until the gas becomes 
hot enough to dissociate H2, leading to the onset of the second 
collapse. 

3.4.1 . The barotropic equation of state approximation 

As mentioned earlier, the easiest way to describe the ther- 
mal evolution of the gas without solving the radiative transfer 
equation is to adopt a barotropic EOS that reproduces both the 
isothermal and adiabatic limiting regimes as a function of the 
density. We use the following barotropic EOS: 



P 
P 



2/3 



(35) 



where p^d is the critical density at which the gas becomes adi- 
abatic. This critical density is obtained by more accurate cal- 
culations and depends on the opacities, the composition and 
the geometry of the molecular gas. At low densities, p <k pad, 
Cs ~ C50 = 0.19 km s"' (for T ~ 10 K), the molecular 
gas is able to radiate freely by thermally coupling to the dust 
and remains isothermal at 10 K. At high densities, p > pad, 
we assume that the cooling due to radiative transfer is trapped 
by the dust opacity. Therefore, P oc p^^^, which corresponds 



In RHD calculations, the radiative transfer equation should 
be formally integrated over 6 dimensions at each time step. This 
process is too computationally demanding for multidimensional 
numerical analysis. Using angular moments of the transfer equa- 
tion is thus very useful, by allowing a large reduction of the com- 
putational cost. However, each evolution equation of a moment 
of the transfer equation involves the next higher order moment of 
the intensity. Consequently, as for the kinetic theory of gases, the 
system must be closed by using an ad hoc relation that gives the 
highest moment as a function of the lower order moments. For 
radiation transport, the closure theory is usually limited to the 
two first moments of the transfer equation. A closure relation for 
the system is needed and this relation is of prime importance. 
Many possible choices for the closure relation exist. In the fol- 
lowing, we present two models based on more or less accurate 
closure relations, that we use in this work: the flux limited diffu- 
sion (FLD) approximation and the Ml model. 

3.4.3. Tine Flux Limited Diffusion approximation 

The diffusion approximation is the most widely used mo- 
ment model of radiation transport. The diffusion limit is valid 
when the photon mean free path is small compared with other 
length scales in the system. On the contrary, the approximation 
is no longer accurate in the transport regime. In the diffusion 
limit, photons diffuse through the material in a random walk. 
Readers can find an accurate derivation of the diffusion limit in 
iMihalas & Mihalas' /l984). §80. 

In the diffusion limit, the radiative energy and radiative 
flux Fr are simply related by 



F,- = V£r, 

3a:r 



(36) 



where /c-r is the Rosseland mean opacity. The radiative flux is ex- 
pressed directly as a function of the radiative energy and is pro- 
portional and collinear to the radiative energy gradient. Equation 
(|36] | has no upper limit, but for optically thin flows, the ef- 
fective propagation speed of the radiation must be limited to c 
(Fr < cEj). We thus have to limit the propagation speed of the 
radiation by means of a flux limiter. Equation (|36] | is then ex- 
pressed as 



cX 

Fr = V£r, 



(37) 



where A = A{R) is the flux limiter. 

In this study, we retain the flux limiter derived by iMinerbol 

(Il978h 



A = 



2/(3 + V9TT2F) if Q<R< 3/2, 
{\+R + Vl + 2R)-^ if 3/2 <R< 00, 



(38) 



with R - |V£'r|/(/fR£'r). The flux limiter has the property that 
A — > 1 / 3 in optically thick regions and /I — > l/Rin optically thin 
regions. 



8 



B. Commerfon et al.: Physical and radiative properties of the first core accretion shock 



In the FLD approximation, an unique diffusion-type equation 
on the radiative energy is thus obtained 



- V ■ (— VEi-l = kp{47tB - cE,). 



(39) 



3.4.4. M1 model 



In the Ml model, the radiation transport is described by 
the first two moments of the radiative transfer equation (ra- 
dia tive energy and flux ) . We u se a closure relation introduced 
by iDubroca & Feugeas' ('1999"). The Ml method is used in the 
RHD code HERACLES (Gonzalez et al. 2007). Based on a mini- 
mum entropy principle, this method is able to account for large 
anisotropy of the radiation as well as for the correct diffusion 
limit. The main advantage of the Ml system is that the underly- 
ing photon distribution function is not isotropic, but has a pref- 
erential direction of propagation. The Ml model also allows to 
explicitly get rid of the ad-hoc limitation of the flux in the trans- 
port regime. 

Let us consider the evolution equations of the zeroth and first 
moments of the specific intensity in the laboratory frame 



d,E, + V-¥, = k^AkB-cE,) 



(40) 



where /cp is the Planck mean opacity. Note that we use the 
Rosseland mean in the first moment equation in order to yield 
the diffusion limit. As a closure relation, the radiative pressure 
is often expressed as Pr = D/Sr, where D is the Eddington ten- 
sor. Assuming that the direction of the radiative flux is an axis of 
symmetr y of the local spec ific intensity, the Eddington tensor is 
given by dLevermord 1 9 84l) 



1 



n igin. 



(41) 



where x is the Eddington factor, I the identity matrix and n = y 
a unit vector aligned with the radiative flux. In the Ml approx- 
imation, the Eddingt on factor y is then obtaine d by minimizing 
the radiative entropy ("Pubr oca & Feugeasll999l) . The Eddington 
factor is then found to be 



3+4/2 



5 -H 2 V4 - 3/2 



(42) 



It is then trivial to recover the two asymptotic regimes of 
radiative transfer. When / — » 0, ;t' = 1/3 and Pi = 1 /3£'il which 
corresponds to the diffusion limit with an isotropic radiative 
pressure. On the other hand, if / = I, X - 1 P, = n ® nE^, 
which corresponds to the free-streaming limit. Between these 
two limits, the closure relation ensures that energy remains 
positive and that the flux is limited (Fj < cE^). 



3.4.5. Systems of RHD equations 

- In the barotropic EOS approximation, the thermal behavior 
of the gas is determined by the choice of the EOS. The en- 
ergy equation becomes superfluous and the Euler equations 
reduce to : 



f 5,p + V [pu] = 

\d,pu + 'Vlpn®n + PI] = -pV<D, 

where O the gravitational potential. 



(43) 



In the FLD approximation, the system of RHD equations 
corresponds to the Euler equations plus the equation on the 
radiative energy 

d,p + V[/)u] =0 
d,pu + V [pu ® u + fl] = -pVO + AVEy 
d,E + V [u (£■ + P)] = -pu ■ V<D + AVE, ■ u 

-KpiAnB - cE,) (44) 

d,E, + V [uE,] = -P, : Vu + V ■ (^V£,) 

+kp(4kB - cE,). 

This system is closed by the perfect gas relation and the 
simplest FLD approximation for the radiative pressure, 

Pr = AE,. 

In the Ml model, we consider the equations governing the 
evolution of an inviscid, radiating fluid 

d,p + V[/5u] =0 
<9,pu + V [pu ® u + PI] - -pVO -I- krFJc 
d,E + V [u (£ + P)] = -pu ■ VO + ^rF,/c ■ u 

-Kp{AnB - cE,) 
d,E, + V [uE,] + V ■ F,. = -Pr : Vu + Kp{4nB - cE,) 
d,¥, + V [uF,.] -H c^V ■ Pr = - (Fr ■ V) u - /fRcF,-. 

where p is the material density, u is the velocity, P the 
isotropic thermal pressure, O the gravitational potential and 
E the fluid total energy E - pe + \ /2pu^. This system is 
closed by the perfect gas relation and the Ml relation (l42t . 

3.5. Initial and boundary conditions 



We use the same model as lMasunaga et al.l d 19981) . i.e. a uni- 
form density sphere of mass Mq = 1 M©, temperature To = 10 
K (Cjo = 0.187 km s^') and radius Rq = lO'* AU. This initial 
setup corresponds to a ratio a of thermal to gravitational ener- 
gies of a = 5RokBTo/i2GMofimu) ~ 0.97 and to a free-fall time 
tff ~ 1.77 X 10^' yr. Boundary conditions are very simple: for 
hydrodynamics, we impose a constant thermal pressure equal to 
the initial pressure (other quantities are free) and for the radia- 
tion field, we impose a vanishing gradient on radiative tempera- 
ture. Calculations have been performed using a Lagrangean grid 
containing 4500 cells. 

3.6. The opacities 

For the Ml model a nd the FLD approxirn ation, we use the 
set of opacities give n by Semenov et alj ( 20031) for low tempera- 



ture (< 1000 K) and'FergusonetalJ ([2005^ for high temperature 
(> 1000 K), that we compute as a function of the gas temperature 
and density. In Fig.|5] the table for the Rosseland opacity is plot- 
ted as a function of temperature and density. In ISemenov et al.l 
(l2003h . the dependence of the evaporation temperatures of ice, 
silicates and iron on gas density are taken into account. In this 
work, we use spherical composite aggregate particles for the 
grain structure and topology and a normal iron content in the 
silicates, Fe/(Fe + Mg)=0.3. 

4. 1D numerical calculations of the first collapse 

As we mentioned before, all calculations presented here have 
been performed with 4500 cells, distributed according to a loga- 
rithmic scale in mass in our initial setup. These mass and spatial 
resolutions are sufficient to resolve quantitatively the accretion 
shock energy budget (see appendix lAll. 
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Fig. 5. Rosseland opacity made of a mix of Semenov et al. 
( l2003h model at low temperature and Ferguson et al.. (,2005,) 
model at high temperature. The Rosseland mean opacity is plot- 
ted as a function of temperature (x-axis) and R 
= ^/10^ (y-axis) using logarithm scales. 



p/(Tl), with 



O 




log(r) 



Fig. 6. Sketch of the estimate of pad used in the barotropic EOS, 
according to the density-temperature distribution obtained with 
the Ml model. 



Calculations with the barotropic EOS have been performed 
using a critical density of p^d - 3.7 x lO"'-' g cm"-' . We derived 
the value of pad from the intersection of extrapolated lines of 
the isotherm and the adiabat in the log(p)-log(r) plane, obtained 
with calculations using the Ml model (see Fig.|6l). 

4.1. Results during the first collapse. Formation of the first 
core. 

Table[T]gives results of barotropic EOS (Baro), Ml, and FLD 
calculations, when the central density reaches Pc = 1 x 10"'° 
g cm"^. Note that for all calculations, the dynamical times are 
very close to ~ 0.189 Myr ~ 1.07 tff. We define the first core 
radius as the radius at which the infall velocity is the largest (i.e., 
the position of the accretion shock). The accretion luminosity 
is estimated at the first core border, according to equation (l34l i. 
The mass accretion rate M and the accretion parameter ^ are 
evaluated at Rf^. Table 1 also displays the central temperature 
Tc, the central specific entropy s^, and the first core temperature 
at the border Tfc . 

As se en from table [H the fir st core radius and mass agree 
well with iMasuna^a et alJ (Il998h results. The results from the 



Ml or FLD calculations are very similar. The central tempera- 
ture with the FLD is slightly higher (by ~ 2%) than the one in 
the Ml calculations, which confirms the fact that the diffusion 
approximation tends to slightly overestimate the cooling. The 
central entropy obtained with the Ml and FLD models is lower 
than the one obtained with the barotropic EOS, and this differ- 
ence increases with time. Indeed, taking into account radiative 
transfer allows the energy to be radiated away during the col- 
lapse, a process which is not properly handled with a barotropic 
EOS. This leads to a lower entropy level of the first core with 
the FLD or the Ml models. The temperature at the border of the 
first core is higher in the Ml and FLD cases, since the photons 
escaping the accretion shock heat up the infalling material. On 
the other hand, the mass accretion rate, the mass, the radius and 
the accretion parameter ^ of the first core are in good agreement 
between the three models. Note that the values of ^ obtained in 
our calculations are closer to the Larson-Penston solution than 
to Shu's expansion wave solution. 

Figure |2] shows the evolution of the central temperature and 
the central entropy as a function of the central density. From Fig. 
13a), we notice the perfect bimodal thermal behavior of the gas, 
from isothermal to adiabatic. The critical density at which the 
gas begins to heat in both Ml and FLD calculations is the same. 
The difference with the barotropic case stems from our choice 
of the barotropic EOS and critical density. For FLD and Ml, we 
find a slope of yeff - 1 ~ 0.61 < 2/3 at high density. This means 
that the first core is not completely adiabatic, but experiences 
some heat loss, yielding a significant cooling. 

In Fig.|3b), the central entropy is plotted as a function of the 
central density. At high density, we recover the adiabatic regime, 
and the Ml and FLD calculations settle at the same entropy level 
in the center, that is lower than the one reached by the barotropic 
model. In the barotropic case, the entropy is determined by the 
EOS. The value of the "adiabat" at which the gas settles is 



i oc In — 



For the barotropic EOS, we have 



i oc In 



„-2/3 2 
P ''sO 



1 + 



Pad 



2/3\ 



oc In 



2/3 
Pad 



P 



,2/3 



(46) 



(47) 



At high density, s — > ^n(c^o/p^[i ^- limiting value for the 
entropy is then ~ 2.03 x 10^ erg K"' g"' for the barotropic 
EOS, which is higher than the value obtained with FLD and 
Ml. Moreover, as the slope of the thermal profile (fig. |3a)) is 
not exactly equal to y = 5/3 in the Ml and FLD calculations, 
as mentioned above, the gas tends to cool and to decrease its 
entropy level at a rate ~ 10"^- lO""* erg K"' g"' s"' at the center. 
This is one of the immediate and most important consequences 
of correctly taking into account radiative transfer in the collapse: 
the accretion shock becomes a real radiative shock and radiation 
is transported outward in the core. 

Figure|8]shows the central temperature, density, entropy and 
optical depth evolution during the collapse. Variations of all vari- 
ables are quite similar for temperature lower than T~ 100 K. 
At 100 K, there is a first discontinuity in the opacity, due to 
the destruction of icy grains. This discontinuity moves to higher 
temperature as density increases, so that density affects to some 
extent the opacity. However, these dififerences are really small. 
From Fig. [He), we clearly see that the entropy level remains 
constant with the barotropic model, whereas the first core keeps 
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Model 






M 






Tk 








(AU) 


(Mo) 


(Mo/yr) 


(Lo) 


(K) 


(K) 


(ergK-' g-') 




Baro 


6.8 


2.35 X 10^' 


4.2 X 10^' 


0.021 


551 


15 


2.03 X 10" 


27 


FLD 


8.1 


2.33 X IQ-^ 


4.1 X 10^^ 


0.019 


411 


58 


2.02 X 10' 


26 


Ml 


7.8 


2.33 X 10^^ 


4.1 X IQ-^ 


0.018 


419 


57 


2.02 X 10' 


26 



Table 1. Summary of first core properties for pe = 1 X 10 g cm 



3.0 L J 2.6x10' 




-20 -18 -16 -14 -12 -10 -20 -18 -16 -14 -12 -10 

log(Po) [g-cm"^] log(pe) [g.cm"^] 

Fig. 7. Evolution of the central temperature ( left) and the central entropy (right) as a function of the central density during the first 
collapse for the Ml (solid line), FLD (dashed red line) and barotropic (dashed-dotted line). 



cooling to lower entropy levels with the Ml and FLD calcula- 
tions. This is an important result, since, as mentioned previously, 
the entropy level of a low mass protostar determines its mass- 
radius relationship, and thus its pre-main sequence subsequent 
evolution. According to the first principle of thermodynamics, 
the global entropy loss of the first core can be estimated as 



[dmj, Uf/,„ \dt)^ \dt 



(48) 



where e is the specific internal energy and v = 1/p is the specific 
volume. If the first core was perfectly adiabatic, the compres- 
sional work would be entirely converted into internal energy, and 
the radiative loss would be zero. As shown in Fig |9f, the lumi- 
nosity increases significantly within the first core in the FLD and 
Ml calculations. This increase corresponds to the radiative loss, 
which amounts to ~ 4x 10"^ Lq. The corresponding entropy loss 
of the first core is then 



(49) 



A.Sfc 1 AL 
~t TAlZ' 

At = 1 X 10-1° g cm-^ we get Mfc ~ 2.3 x lO^^ Mq. The 
entropy loss is thus 



Asfc 
At 



3.3 X 10- 



erg K 



(50) 



For a typical first core temperature of a few 100 K, the entropy 
loss rate is ~ 10"^ erg K"' g"' yr"', which is consistent with the 
value found at the center The first core lasting about a few hun- 
dred years, the total entropy loss thus represents a few percents 
of the core's initial entropy content, a nd this entropy l oss in - 
creases with time, as shown in Fig. 1 of iMasunaea et al.l (Il998h . 
Such an entropy loss cannot be handled with a barotropic EOS. 

In Fig. [8jd), we see that the first core quickly becomes op- 
tically thick once it begins to heat up. The optical depth at the 
center is so high that observations are not able to catch the cen- 
tral evolution of the cores at this stage of the evolution. 



4.2. Mechanical and tiiermal profile of the first prestellar 
core. 

Figure |9] displays the radial profiles of the density, gas tem- 
perature, velocity, specific entropy, optical depth, luminosity, ra- 
diative flux and integrated mass, once the central density reaches 
Pc = 1 X lO""' g cm"^, for calculations with Ml, FLD and the 
barotropic EOS. As mentioned before, the first core radius is 
smaller with the barotropic model, whereas both Ml and FLD 
results are similar For the latter models, we find only small dif- 
ferences around r ~ 1, i.e. at the transition between the optically 
thick and thin regimes. Since the Ml model naturally recovers 
the diffusion and transport regimes and since the FLD model 
is defined as to recover these limits as well, it is natural that 
we get similar results when either the transport or the diffusion 
regimes are well established. Although the accretion shock is lo- 
cated within the transition region around r ~ 10, the aforemen- 
tioned small differences do not affect the first core properties. 
The entropy jump at the accretion shock is much higher with 
Ml and FLD than with the barotropic approximation. We also 
see from the temperature and entropy profiles that the barotropic 
EOS cannot handle correctly the transition from the isothermal 
to the adiabatic regime, as pointed out earlier. The radiative pre- 
cursor in front of the shock is not reproduced with the barotropic 
EOS. In that case, the gas becomes rapidly isothermal after the 
shock, whereas it is heated by photons escaping from the shock 
in the FLD and Ml models. As a consequence, the core tem- 
perature is higher in the barotropic case. The differences in the 
behavior of the radiative flux at small radius come from the dif- 
ferences in temperature between the Ml and FLD models. 

The emergent luminosity, on the other hand, is the same for 
both the Ml and FLD models. The luminosity jump, ~ 0.013 
Lq is consistent with the accretion luminosity estimated from 
equation (|34] |. 0.014 L©. This means that all the infalling kinetic 
energy is radiated away at the shock boundary, i.e. that the shock 
is supercritical at the formation of the first core. 
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Fig. 8. Evolution of the central temperature, density, entropy and optical depth with time for the Ml (solid line), FLD (dashed red 
line) and barotropic (dashed-dotted line) models. 



If we apply the criterion derived in equations ([20} and (|211 
to the upstream quantities estimated in Fig.|9] i.e. pi ~ 7 x 10"'"^ 
g cm""*. Ml ~ 2 X let'em s"', and Ti ~ 57 K, the corresponding 
critical temperatures are T„ ~ 23 K with ( |20| ) and ~ 57 K 
with i2T[ . This confirms the fact that the accretion shock on the 
first core is supercritical. However, since accretion takes place 
in the transition region between optically thin and thick regimes 
and since most of the upstream region is optically thin (see lu- 
minosity and optical depth profiles in Fig. |9]l, the most relevant 
model is the one we developed for a supercritical shock with a 
upstream optically thin material. 



4.3. Comparison with an analytical model 

In Sect. 2, we have developed a semi-analytical solver 
that can be applied to the accretion radiative shock on the 
first core and that covers three cases: the case of sub- or 
super-critical shocks in an optically thick material and the 
case of a supercritical shock in an optically thin material. In 
the following, we develop a simple model for the protostellar 
collapse, where the upstream quantities are estimated under 
some basic assumptions, and we compare this model with 
our numerical results. We need to estimate the density, the 
velocity and the temperature in the preshock region in the 
context of the first core formation. We can easily get the 
density (and the velocity) from the Larson-Penston and Shu 
self-similar solutions. The tricky part is to estimate the tem- 
perature before the shock. However, as mentioned in ^2.41 the 
preshock temperature is determined by the upstream velocity. 
Assuming that the accretion shock is supercritical and that 
the upstream material is optically thin, it is then easy to get 



this te mperature and to recover all fluid variables. lOmukail 
(2007) proposes an alternative model that fully describes the 
first core characteristics but does not use jump relations for a 
radiating fluid. In our model, we consider the characteristics at 
the first core border and the jump conditions for a radiating fluid. 



We approximate the upstream velocity by the free-fall veloc- 



ity 



Ml.ff - 



2GMk 



(51) 



The preshock density is given assuming a profile oc r ^ in the 
free-falling envelope 



Pi 



-2 



(52) 



where Cj = (P/py^^ is the isothermal sound speed (at 10 K). The 
temperature is estimated by assuming a supercritical shock, i.e. 
all the upstream kinetic energy being radiated away. Moreover, 
our calculations using the FLD or Ml models yield M\ ~ 2 and 
Fig. |2] shows that X ~ I. The shock temperature thus reads (c.f. 
equationl22l) 



piu 



I.fF 



(53) 



We take the first core properti es obtained in our num eri- 
cal calculations, similar to those in Masun aga et al.l (11998 *). i.e. 
Mfc ~ 2.3 X 10"^ Mq and Rf^ ~ 8 AU. If we apply equations 
(ISTT l. (l52T i and (|53]l to these quantities, we get: pi - 5.8 x 10"'"* 
g cm-^ Ml ~ 2.3 X 10^ cm s"' and Ti ~ 50 K. 
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Fig. 9. Radial profiles of various properties during the collapse of a 1 dense clump for a core central density pc - I x 10"'" g 
cm"^, for the Ml (solid line), FLD (dashed red line) and barotropic (dashed-dotted line) models. From top left to bottom right: (a) 
density, (b) gas temperature, (c) velocity, (d) entropy, (e) optical depth, (/) luminosity, (g) radiative flux and (h) integrated mass. 



We now apply these preschock quantities to the jump prop- 
erties derived in Sect. 12.2.21 for a supercritical shock with an 
optically thin upstream material. According to eq. [14] we get 
- Pi/pi ~ 50, which gives p2 ~ 2.9 x 10"'^ g cm"-', and 
M2 = 4.6 X 10^ cm s"' . Using equation ( fTSl l, we get X ~ 0.99: aU 
the infalling kinetic energy is radiated away at the shock. 



We can compare these analytical results with the numerical 
ones, given in Fig. |9] and in Table [T] pi ~ 7 x 10"'"* g cm"^. 
Mi ~ 2 X 10^ cm s"' and Ti = 57 K. Thus, the analytical 
estimates for the preshock quantities agree quite well with 
the numerical values. We also read from the figure and table 
P2 ~ 2.15 X 10"'^ g cm"^, M2 ~ 10"* cm s"' ui, again very 
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similar to our analytical estimates. These comparisons validate 
our simple analytical model to infer the properties of the first 
Larson's core. 



5. Summary and perspectives 

In this paper, we have conducted RHD calculations aimed 
at describing the physical and radiative properties of radiative 
shocks. We have applied these calculations to the formation of 
the first Larson's core, during the early phases of protostellar col- 
lapse. We have also derived simple analytical models and com- 
pared the results with the ones obtained in the simulations. The 
main results of our study can be summarized as follows : 

1. The properties of the first collapse and the characteristics of 
the first core in our calculations are i n good agreement with 
the ones obtained bv lMasunaga et al]([l998i) . We find that the 
first core has a typical radius of ~ 8 AU and a mass ~ 2 x 1 0"^ 
Mq. This sets up the initial conditions for the second collapse 
and the formation of the second Larson's core. 

2. We show that, at the first core stage, the accretion shock is 
a supercritical radiative shock, at which all the infalling ki- 
netic energy is radiated away, and that a barotropic EOS can- 
not reproduce the coiTect jump conditions at the shock. The 
FLD and Ml calculations show that there is a substantial en- 
tropy loss during the formation of the first core, due to the 
radiative loss. A barotropic EOS cannot handle correctly this 
cooling mechanism. In consequence, the first core's entropy 
content obtained with a barotropic approximation is over- 
estimated, compared with the calculations which solve the 
radiative transfer. Such a cooling effect can have a strong 
impact on the core fragmentation process and, eventually on 
the initial properties of the future protostar (Commergon et 
al., in prep). 

3. We confirm that, when radiative cooling is properly taken 
into account, the transition from an isothermal to an adiabatic 
phase during the first collapse and the formation of the first 
Larson's core does not necessarily correspond to an optical 
depth of unity, as shown initially by Masunaga & Inutsuka 
(1999). 

4. We develop a simple analytical model for supercritical 
shocks within an optically thin medium, which reproduces 
well the jump conditions obtained with the numerical calcu- 
lations. We show that the compression ratio in such a kind 
of shock can become very high (r ~ 50). We plan in the 
future to keep exploring this issue with a frequency depen- 
dent radiative transfer model. Indeed, strong shocks on (mas- 
sive) protostars are known to be optically thick for hard pho- 
tons, while optically thin for UV radiation (e.g. IStahler et al.l 
119801) . 

5. We show that, at least in ID spherical calculations, the flux 
limited diffusion approximation is appropriate to study the 
earliest stages of star formation, as it gives very similar re- 
sults as the more complete Ml model for radiative transfer. 
Note, however, that our ID spherical geometry code cannot 
account for multi dimensional effects like the anisotropy of 
the radiation field. 

This study confirms the necessity to solve the complete RHD 
equations, i.e. to correctly take into account the coupling be- 
tween radiation and hydrodynamics, when addressing strong 
shock conditions, as occurring during the formation of the first 
Larson's core in the context of star formation. Such a complete. 



consistent treatment of radiation and hydrodynamics is neces- 
sary to correctly handle the cooling properties of the accreting 
gas and thus to obtain the correct mechanical and thermal prop- 
erties of the first core. This, in turn, sets up the initial conditions 
for the second collapse and the formation of the second core. 
Indeed, since the first core has a short lifetime (a few hundred 
years), it is important to pursue these calculations during the 
second collapse, including H2 dissociation. Such work is under 
progress. 
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Appendix A: Effect of the numerical resolution on 
the energy budget through the shock. Case of a 
0.01 Mq dense core. 

In contrast to the hydrodynamical case, the structure of a ra- 
diative shock can extend over large distances, depending on the 
optical properties of the material. For an optically thin material, 
the photon mean free path is large, so the shock structure is very 
extended compared to the viscous scale. In this work (see Sect. 
O, we present numerical calculations of dense core collapse, us- 
ing a fix resolution in mass, i.e. the mesh is not refined in the 
large gradient zones. Although this Lagrangean description is 
well suited for the hydrodynamical shocks, thanks to the artifi- 
cial viscosity scheme, it may encounter difficulties in the case 
of radiation-hydrodynamical flows, in particular in the optically 
thin region (upstream region, outside the first core). 

In this appendix, we present the results of the collapse of a 
0.01 Mq dense core, using the same initial ratio of thermal en- 
ergy over gravitational energy as in Sec|3](Q' ~ 0.97). To investi- 
gate the effect of the numerical resolution, we performed calcu- 
lations with 4500 cells and 18000 cells, using the FLD model. 

Figure IATI shows the density, gas temperature, velocity, en- 
tropy, optical depth and luminosity radial profiles for the 2 cal- 
culations at a central density pc ~ 1.6 x 10 " g cm"^. Although 
there are some significative differences in the radiative precursor 
region (i.e. the transition region between optically thin and thick 
regions, where 2 < t < 0.5) and in the estimate of the first core 
radius (~ 10%), the entropy, density, velocity and luminosity 
jumps are about the same. In both calculations, the shock is su- 
percritical and the amount of energy radiated away is about the 
same (L = 1.5 x 10"^ with 4500 cells, and L = 1.44 x 10"- with 
18000 cells). This means that the overall properties of the first 
accretion shock, including its global energy budget remain cor- 
rectly calculated even at low resolution. However, using 18000 
cells, we see that the spike in the gas temperature is resolved and 
that the radiative precursor length is much smaller. On the other 
hand, the central entropy within the first core is the same in both 
cases, indicating that the cooling of the first core by radiation is 
not affected by the lack of resolution in the radiative shock. 

Figure IA.2I shows the normalized energy balance at pc - 
1 .8 X 1 0" " g cm"-' for the calculations run with 1 8000 cells (left) 
and 4500 cells (right). The figures display the rate of change of 
potential energy dEp/dt, kinetic energy dE^/dt, internal energy 
dUJdt, total energy d(E\^ + E-p + U\)ldt, and the work done by 
thermal pressure and radiative flux (L^d + A-nr^Pu). The total 
energy equation reads: 

-[U, + E^ + Ep + E,) + — [Anr\u{p + A) + F,]] = 0. (A. 1) 

First, we see form fig. lA. 21 that the radiative pressure exerts a 
negligible work compared to the thermal pressure. Comparing 
the energy balance of the two calculations, we see that it is glob- 
ally the same, which confirms that the calculations done with 
4500 cells get the correct features of the first core accretion 
shock, even though the radiative structure of this shock is not 
resolved. 
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Fig.A.l. Radial profiles of various 1st core properties during the collapse of a 0.01 clump for a core central density pc - 
1 .8 X 10"" g cm"'', for calculations done with 4500 cells (dotted red line) and 18 000 cells (solid black line). From top left to bottom 
right: {a) density, {b) gas temperature, (c) entropy, (d) velocity, (e) optical depth, (/) luminosity. 



d{E^ + + U)/dt 

d(E,)/dt 
d(EJ/dt 
d(U)/dt 




-a 



-4 -3 
log(M) [M,„„] 



d(E^ + + U)/dt 
L„a + 47rr=PjU 
d(E,)/dt 
d(Ej/dt 
d(U)/dt 




-4 -3 
log(M) [M J 



Fig. A.2. Normalized energy balance as a function of the mass for the calculations done with 1 8000 cells (left) and 4500 cells (right), 
when the central density reaches pc - l.^ x 10"" g cm"^. Scales are logarithmic and normalized to the rate of change of the total 
energy d(E]^ + Ep + Ui)/dt. 



